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NASA Grant Number NSG-1248 



The NASA Technical Officer for this grant is Mr. Manuel Salas, 
NASA Langley Research Center, Hampton, Virginia. 


In the first semester of the one-year grant, the following tasks have 
been performed: 

(1) Formulation of a generalized method of conformal mapping to 
generate computational grids about cross-sections of actual airplanes, 
including the case of wing cross-sections detached from fuselage cross- 
sections. The method is entirely analytical and it consists of a repeated 
application of the Karman-Trefftz mapping function. Techniques to make 
the application of the method automatic have also been developed. The 
pertinent FORTRAN codes are ready to be used. The analysis has been 
reported in a paper accepted at the Symposium on Numerical Laboratory 
Computer Methods in Fluid Mechanics, sponsored by the ASME , in New 
York, November 1976, and which is reproduced in full in the following pages. 

(2) Formulation of the geometry of a sample case, a simplified arrow-wing 
airplane. A FORTRAN code to generate such a geometry at any number of 
specified cross-sections and with a variety of basic parameters has been pro- 
vided to NASA on February, 1976. 

(3) Formulation of the gasdynamical analysis of a three-dimensional, super- 
sonic, inviscid, steady, shockless flow past an arbitrary airframe, using 
the computational grids generated by (1). The analysis, which includes special 
treatments for body points and bow-shock points, relies on the equations of 
motion written in terms of logarithm of pressure, entropy and two angles re- 
prising the velocity vector, A FORTRAN code including the sample geometry 
of (2) has been prepared and is almost totally debugged. To date, the flow 
analysis is considered reliable wherever the cross-section is elliptic in shape. 
Further tests are being conducted to determine reliability in more complicated 
cases. A detailed report on the analysis is in the process of being written and 
it will be submitted to NASA shortly. Results of the calculations will be sub- 
mitted for acceptation at the winter meeting of AIAA, January 1977. 


CONFORMAL MAPPINGS FOR COMPUTATIONS OF 
STEADY. THREE-DIMENSIONAL. SUPERSONIC FLOWS 


1. Introduction 

Steady, supersonic, three-dimensional, inviscid flows can be 
evaluated numerically by marching techniques, which amount to integrating 
the equations of motion in successive steps along a fixed axis lying in the 
general direction of the velocity vectors. In the numerical analysis of the 
flow about a vehicle, we will stipulate that our basic frame of reference is 
a Carte sian- (x, y, t) frame, whose t-axis lies along the body and whose 
(y, t) plane is the plane of symmetry of the body. The symbol t is used 
here in lieu of the conventional z to avoid confusion with z as a complex 
variable; the choice of t is suggested by the time-like role which t plays 
in the numerical procedure. In fact, at each computational step a set of 
initial values is known at a certain number of points on a plane, normal to 
the t-axis and a set of final values is to be determined at corresponding 
points on another plane, also normal to the t-axis, downstream of the ini- 
tial plane; At is the distance between the two planes. 

For the sake of accuracy, the following conditions should be satis- 
fied; 

1) Computing at evenly spaced grid points; 

2) Having a set of grid points exactly on the body and on the outer 
boundary (the bow shock); 

3) Making the grid in each t = constant plane as close as possible 
to an orthogonal grid, particularly in the vicinity of the body. 

The Cartesian frame being obviously unable to satisfy the second 
condition, it is necessary to recast the equations of motion with reference 
to another set of independent variables, X, Y and T, such that the X=0 


plane represents the body surface and the X=1 plane represents the bow 
shock wave; in addition, to satisfy the third condition, the two families 
of X= constant and Y = constant lines in any plane normal to the t-axis 
shovild be as close as possible to an orthogonal net, at least in the vicinity 
of the body. 

If the geometry of the rigid boundaries is simple, the mapping can 
be obtained in a straightforward manner. For example, to study the flow 
past a nearly circular cone (Fig. la), the equations of motion can be writ- 
ten in a cylindrical (p, 9, t) frame, with 

X s p cos 6 

(1) . y = p sin 0 

t = T 

Let p = b(0,T) be the equation of the body wall and p = c(6,t) be the equa- 
tion of the bow shock. Then, the transformation 



(2) Y = 0 

T = T 

maps the shock layer onto a parallelepiped, a cross section of which is 

s|e 

the rectangle shown in Fig. lb. In it, the image of the body is a straight 
segment, AC (X=0), the image of the bow shock is a straight segment, 

DF (X=l) and AD and CF are symmetry lines (y s ± n/2). Lines defined 
by evenly spaced values of X and Y are shown in both the computational 
and the physical plane. The equations of motion are, once more, easily 

Stretchings in either one of the X- and Y- directions can be used, if 
necessary; in which cases Equations (2) are replaced by more compli- 
cated expressions. We are not considering such possibilities in the 
present paper, since our attention is focussed on the part of the mapping 
problem, so far expressed by (1), 
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recast in terms of the computational variables, X, Y and T. 

Note that , if the cross sections of the body and the bow shock were 
concentric circles, the grid in the physical plane would be orthogonal, but 
this is generally not the case. The lack of perfect orthogonality, however, 
does not impair the accuracy of the results to within reasonable limits. 

To maintain accuracy in the computation of the flow about an 
elliptic cone, with an axis ratio, say, greater than 5, a more sophisticated 
mapping is needed. In fact, if we proceeded as in the preceding case, the 
grid in the physical space (shown in Fig. 2a) would be grossly inadequate 
in the vicinity of the body. Elliptic coordinates should be adopted in each 
cross section^, resulting in the grid shown in Fig. 2b. 

The construction of such a grid depends on a simple conformal 
mapping. Let 

(3) z = X + iy 

and 

(4) C = 0e‘® 
be two complex variables, and let 

(5) z=j(C+^) 

where h is a real number. If A and B are the semi-axes of the ellipse 
in Fig. 2b, and 

(6) h = (A2 _ b2)^ 

the ellipse is mapped onto a circle in the C-pla>ie, having its center at the 
origin and a radius, p^, defined by 



The portion of the z-plane exterior to the ellipse is mapped onto the 
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portion of the C-plane exterior to the circle. Obviously, the image of the 
shock in the C-plane is another closed line (all the more similar to the 
shock itself the farther it is from the body since, at infinity, z tends to j C). 
Let it be called the p = c line. Then, one can proceed as before, using 
( 2 ). 

In what follows, I intend to generalize the concept above, by showing: 

1) that simple rules can be found to map the portion of a plane exte- 
rior to a given rigid body (or bodies) conformally onto the portion of a plane 
exterior to a nearly circtilar contour, so that the final (non conformal) 
mapping described by (2) and illustrated in Fig. 1 can always be applied; 

2) that such mappings can be obtained by repeated applications of 
a single, and rather elementary, analytic function; 

3) consequently, that the application can be automatized, without 
requiring any guesswork by the user; 

4) that all derivatives required for the numerical analysis of the 
flow can also be evaluated analytically and automatically. 

A glance at Fig. 3, where two typical cross sections of an airplane 
are shown, gives an idea of how much bigger the challenge is in dealing 
with practical applications than in the idealized case of the elliptic cone.- 
It also shows how drastically the contour to be mapped changes from one 
cross section to another of the same airplane. 

2. Numerical Mapping Technique 

Obviously, one could solve the conformal mapping problem numer- 
ically on each and every cross sectional plane. The physical shock layer 
could even be conformally mapped directly onto a rectangle, in this way 
eliminating the second, non-conformal step expressed by (2). In principle, 
the problem could be solved by simply adding to the equations of motion 
the two equations: 
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( 8 ) 



which assert that the complex variable Z = X 1 iY is an analytic function 
of z defined by the condition that the physical boundaries must be mapped 
onto the sides of the rectangle, as shown in Fig. 1. The feasibility of 
such an approach has been demonstrated by J. F. Thompson, et al. in a 
recent paper^, which contains a host of seemingly successful applications 
to incompressible flow problems. Any attempt to compare the strictly 
numerical approach to the one described in the present paper would be 
premature. A comparison will requi/e a good amount of additional informa- 
tion; definitive data on computational time are not available; neither do we 
know, at this time, whether the number of grid points needed to solve the 
mapping problem numerically would be compatible with, or would exceed 
the minimum number of grid points needed to perform the flow field analy- 
sis. For a three-dimensional, compressible flow, the equations are much 
more complicated than for a two-dimensional, incompressible flow, even 
if the latter is viscous. The role played by time in the latter and by the 
time -like variable t in the former is different; more 
t-dependent terms appear in the latter than in the former. A possible 
advantage of the present technique may stem from the fact that it defines 
the mapping and all related derivatives analytically—that is, exactly. 
Derivatives produced by the gradual variation in the mapping function along 
the t-axis, as a result of which the physical grid changes whereas the 
computational grid remains unchanged, are also definable analytically. It 
may be worthwhile to mention that full-scale computations of three-dimen- 
sional, steady, supersonic flows past airplanes using conformal mappings 
have been performed already® without lengthening the computational time by 
more than a factor of two. 
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3, Analytical Mapping Technique 

As a first step towards constructing a mapping function, we observe 
that a typical cross section (such as in Fig. 3a) contains a large number of 
corners and edges. Each sharp angle in a contour (and its symmetrical on 
the left half plane) is smoothed out by a mapping of the type: 



where h is the complex coordinate of the corner or edge in the z-plane, 
h is its conjugate and 6 is the inverse of the external angle in multiples 
of TT. In the f -plane, T = ± 1 are the images of the two points, z=h, 
z=-h but the contour is now smooth at C = - 1* All other sharp angles on 
the contour remain unaltered in the Q-pla.ne. The mapping defined by (9) 
is, to within a rcale factor, the classic Karman-Trefftz airfoil -gene rating 
function. Note also that the function defined by (5), known as a Joukowski 
mapping, 's a particular case of (9) with a real h and 6 = 1/2. 

The same procedure can obviously be applied iteratively, each 
time with a different value of h and a different value of 6, until all cor- 
ners and edges have been eliminated. Since (9) is symmetric with respect 
to the imaginary axis, the successive contours remain symmetric in the 
intermediate planes and in the final plane. The final contour is certainly 
smooth and, in all cases of practical interest considered so far, very close 
to a circle. 

4. Additional Remarks 

The greatest difficulty in writing a general purpose code for the 
repeated application of (9) stems from the multi-valuedness of the arc- 
tangent and the need for continuity in the mapping, except across the cut 
running between the branch-points of (9). Note, first, that (9) can be 
decomposed into two parts: 
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In turn, (10) can be computed in the form, 

(12) d = ^ ~ - - exp{i6[arg(z - h) - arg(z + h*^)]! 

z + h 

Now, 

(13) arg(z-h) = arctan arg(z+h*) = arctan = Qp 

if 

(14) h = a + iP 

The angles 6j^ and 6p are shown in Fig. 4 where the heavy line represents 
a possible cut. It is clear that, as z moves in the right half plane, 6p 
varies continuously and it can be correctly defined by standard ATAN 
functions in Fortran. The case for is more difficult. So long as z 
remains to the right ofAB, where Real (z-h) > 0, the Fortran definition of 
6j^ is correct. If z enters the region to the left of AB, below the cut, -n 
must be added to as defined by Fortran; if that happens above the cut, 

TT must be added. In Fig. 5 more extravagant cuts are shown. In addition 
to the regions, similar to the ones discussed above, where ±n should be 
added to 6j^, there are two possible regions, shadowed in the figure, 
where ± 2n should be added, in order to assure continuity in 0j^. In 
coding, the difficulty can be circumvented very easily, by introducing a 
factor, M, equal to 1 for all points above the cut and to -1 for all points 
below the cut. If 6^ is the value of arctan [(y-P)/(x-tt)] as provided by 
Fortran, 
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Real (z-h) > 0, sgnM = sgn[lmag(z-h)] 


, Real (z-h) < 0 

\ + 2utt , Real (z-h) > 0, sgnM 4 sgn[lmag(z-h)] 

The concepts of "above the cut" and "below the cut" are invariant by con- 
formal mapping; therefore, the value of U for a given point with respect to 
a given cut need be defined only once even when a string of intermediate 
mapp^.-gs is used; this can be done, thus, in either one of two planes where 
the dellnition is simple, that is, either on the z-plane or on the final 
C-plano (the plane of the near-circle). In Fig. 6, for example, the cuts 
needed in the successive intermediate mappings which transform the shape 
of Fig. 3a into a near-circle are shown, both in the z-plane and in the 
C -plane. The definition of the factors U in the C -plane is obvious. 

In principle, the intermediate mappings can be e:.;ecuted in any 
order. In practice, the coding turns out to be much simpler (and, in many 
cases, the contour in the C-plane closer to a circle) if the edges and cor- 
ners are eliminated, beginning with the one having the smallest value of 
6 (largest angle) and in the order of increasing 6’s, For example, in the 
case of Fig. 6, if the corner at D were eliminated before eliminating the 
edge at E, the corner at F would disappear simultaneously and a special 
logic should be provided to take care of such a possibility. If the edge at 
E is eliminated first, F and D are far apart from each other in the fol- 
lowing steps and can be treaced as usual. An even more delicate situation 
would arise for a cross section of the type shown in Fig. 7, unless the 
edge at B is eliminated before treating the corners at A and C. Indeed, 
if the corner at A (having an exterior angle equal to n/3) were to be 
eliminated first, part of the shock layer would penetrate into a second 


(15) 


M 


00 
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Riemann sheet (u'le to the fact that 6=3) and should be recalled back 
when treating C. It is clear, instead, that ordering the edges and corners 
by increasing values of 6 results in a simple logical operation. 

To conclude our general considerations, let it be noted that, 
whether we want to proceed from the z-plane to the C*-pl&ne or vice versa, 
the nature of the intermediate mappings is the same. In fact (9) can be 
written as follows: 



and (16) can be used to obtain z from ^ as (9) was used to obtain C from 
z. Therefore, as announced above, the mapping problem can be reduced 
to repeated applications of a function, MAP, a sample coding of which is 
given in Fig. 8, regardless of the complexity of the geometry and of the 
direction in which the mapping is considered(from the physical cross sec> 
tion to the near-circle, or vice ve' a). In the calling sequence, INDEX 
is set equal to 1 if (9) is used to get C from z (in v'hich case POINT 
means z and MAP means Q), equal to -1 if (16) is used to get z from Q 
(in which case POINT means ( and MAP means z); FACTOR is any 
convenient number havx. :he same sign as u. The meaning and use of 
DER will be discussed later on. Note also that H(J, J) and POWER(J) are 
h and 6 respectively as they appear in (9) and (16). Finally, note that in 
the program of Fig. 8 the special cases shown in Fig. 5 are not contem- 
plated. 

5. Examples 

We turn now to examining few examples of mappings, involving 
non-trivial geometries. Note that the input is the geometry alone (which 
obviously includes the location of edges and corners in the physical plane). 
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The program orders the edges and corners in order of increasing ft, and 
finds the location of their images in the C-pl^ne. Then, at regularly 
spaced values of 6 between -n/2 and tt/2, it finds, by trial and error 
(and using the mappings in reverse), which value of b produces a point on 
the original contour in the z-plane; generally, three iterations are suffi- 
cient. In these test cases, a shock is assumed as a circular boundary in 
the z-plane. Once more, by trial and error, all the proper values of c(6) 
are found. Finally, for each 0, and for intermediate values of X, de- 
fined by 

07) = (l^J 

the remaining points on the grid are determined. Note the slight difference 
between (17) and the first of (2), which permits grid lines to concentrate in 
the vicinity of the body. The figures are self-explanatory. In Fig. 9 the 
case of Fig. 7 is solved; in Fig. 10 we see a similar problem (with equal 
values of the corner angles) but with a more clustered geometry; in Fig. 11 
a case similar to the one of Fig. 3a is considered (an airplane with a fin on 
the wing); and, finally, in Fig. 12 a cross section of a fuselage and an 
arrow wing is mapped. Note that in the latter case the near-circle contour 
is composed of alternating parts representing rigid walls and points in the 
free flow. 

6. Coupling of Mappings and Equations of Motion 

In non-dimensional form, let P be the logarithm of pressure, S 
the entropy, V? the temperature, defined by 

(18) ^=exp(;i-^P+is) 

and V the velocity vector. The equations of motion are: 
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V. 7p + Y^. V = 0 

(19) - Vx V X V +^Vp = 0 

V.7S = 0 

A A A 

Let I, J, K be the unit vectors in the direction of the x-, y-, and t-axi**, 

A A 

respectively. Let also i and j be the unit vectors in the z-plane, parallel 
to the 6=:constant line and the p-constant line at each point. Let us denote 
by an index, all vector operators defined in a t=constant plane only, 
and let 

(20) V = w(y + K) 

It is e af.y ro prove that the equations of motion (19) can be recast in the 
form: 

1 7/r) - XX ^.x r- p (7.p-?a ) + 

l.7,S =■ 0 

Scalar forms can be easily obtained from (21) but it must be kept 
in mind that, if non-Cartesian variables such as p and 9 are used in each 
cross-sectional plane, i^e values of x and y at a node of the grid (at con- 
stant p and 6) change from plane to plane; therefore, for any function, 7: 



(22) 

?t = ^p ^ ®t 

For 

example, if 

^ A ^ 


(23) 

V = 0 i + T|y 

and 

• (24) 

dC ^ iuu 
g = — 2. = G e 
® dz 


where C is the complex variable of the near-circle plane, obtained, if 
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necessary, through repeated applications of (9), it follows that 


Let 


x = f V - v:i (’li 


(2A) 

£ = f, +if. =i^ 

(27) 


(28) 

'I' = 't'l + ii'2 = ^ 

(29) 

A, =g2 + pfj. A, =5i+f,, [y 

(30) 

= Ga+pfj, Ba= — Ti + f0 

(31) 

D = ^ [ti(1 - coi) - a CD®] + fa - ilf# 


>--) 


Then, (21) can be written in the form; 

<7^ + BjCTp + BjCT^ + ^ 
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A final change to the variables X, Y and T can now be performed but will 
not be considered in the present paper. As one can see by inspecting (32), 
prior to integrating the equations of motion one must evaluate g, f, cp, and 
i|f. In what follows, we will briefly expose how the evaluation is performed 
in a complicated case. 

Let us assume that J successive mappings are necessary to elimi- 
nate all edges and corners, each partial mapping being of the form (9), 


that is: 


2 ,-^. - 1 




. -+-1 ' ' 


We will also identify the physical z-plane with the Zj -plane. In such a 
plane we consider J+2 points, denoted by h^j (Fig. 13); the first J points 
are representative of edges and corners and will operate as "hinges" in 
the successive mappings; the last two points are the intersection of the 
physical contour with the y-ajcis. The image of h^j in the z^ -plane will 
be called h^^. Note that h^j is the active hinge in the j-th transformation 
and, therefore, 

(34) h. = 1 

The J-th mapping produces a smooth contour, but it is convenient 
to have it as centered as possible about the origin of the There- 

fore, a (J+1) st mapping is used, which is simply a translation: 


C = Z, 


S , S = y (h 


I'^^J+l, J+1 ^J+2, J+1^ 


In conclusion. 




For each partial mapping except (35), let 
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Dt 


~ - Lit “j it ^ i aj+i^ J 


Similarly, 


2 ^ 5i 


5i &t 


J*2.0*l 


Note that Bh. ..i/^t can be computed as a particular case of (44): 
j+i 






Bh. . 


Note also that, if -t=j, — — = 0, a result consistent with (34). 

The values of Bh. ,/Bt (-1=1, J+2) and of B6./Bt are known from the 
vl J 

geometry of the vehicle as a function of t. Therefore (46) permits all the 

values of 9h. . i/^t (j=l through J) to be computed recursively and (44) 
J+ i » J+i 

furnishes all the values of Bz^^j/Bt, including Bzj^, /Bt, in a similar way. 
Then, BC/Bt is obtained from (45) and f is obviously equal to . 

From Fig. 8 it appears that g^, as defined by (37), is computed 
every time either Zj^j(Zj) or z^(Zj^j) is computed, and in either case it 
has the same meaning. The factor, g^, appearing in (46), must be under- 


stood as 


Xj±L 


and is thus obtainable when MAP is applied to determine the values of 

h, from the values of k,.. 
t,j+l Ij 

Finally, to compute ii' one uses (39) and (37) obtaining, after more 
manipulations. 
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(47) 


2zj,i azM 


-2 




ilj 


At this stage, all elements necessary to evi. ■ ate (47) numerically 
are available. 

7, Conclusion 

Sections 4 and 6 contain all the information which is necessary to 
use the conformal mapping concept, as explained in Section 3, for the 
numerical analysis of a three-dimensional, steady, supersonic, inviscid 
flow. Details of such analysis, from the fluid -mechanical viewpoint, and 
a sample calculation of a redime ntary arrow-winged plane will be pre- 
sented in a forthcoming Report. 
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COMPLEX FUNCTION MAP ( iNDEXt J*POINTtFACTOP*DER) 

COMPLEX POINT f HINGE tCI«P*0«D£R«H( 10 • 10) 

C0MM0N/CMAP/PI*PI02*CI »HtPOWEP(10) 

C •••• DER IS DERIVATIVE D (Z (J« 1 ) ) /O ( Z ( J) ) •••• 

HINGE=H(J*J) $ POW=POWE«(J) S IF ( INDEX .EO.-l ) POWsl ,/POirf 
PsPOiNT S Dl=PEAL(HlNGt) S 02=AlMAG (HINGE ) 

IF ( Index. eo,1)p=(point-h Inge )/D i»i, 

XM=REAL(P)-1. S YM=AIMAG(P) $ XP=XM42.S Y2aVM*»2 S RM2*XM**2*Y2 
RP2=XP**2»Y2 5 THP=ATAN(YM/XP) % IF (XM.EO.O. )GO TO 1 
THM=ATAN(YM/XM) $ IF(XM.6T.0.)G0 TO 2 S THMrTHM»SI6N(PI*FACT0R) 
GO TO 2 

1 THM=SIGN(PI02*YM) 

2 D=(Rm2/RP2)«*(.5*P0»^)*CEXP(CI*P0W*(THm-THP) ) 

MAP=(1.40)/(1.-D) $ IF(INOEX.CO.-l)GO TO 4 
DER=D1®P0W*(MAP**2-1.)/( (POINT-HINGE )*(P0INT*CMPLX(D1*-D2) ) ) 
RETURN 

4 MAP=HINGE*(MAP-1.)*D1 

OER=DI •(P0INT**2-1.)/(P0W* (MAP-HINGE )* (MAP4CMPLX (D1 f-D2 ) ) ) 
RETURN $ END 












